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ABSTRACT 

We investigate the gravitational fragmentation of expanding shells in the context of 
the linear thin-shell analysis. We make use of two very different numerical schemes; the 
FLASH Adaptive Mesh Refinement code and a version of the Benz Smoothed Particle 
Hydrodynamics code. We find that the agreement between the two codes is excellent. 
We use our numerical results to test the thin-shell approximation and we find that the 
external pressure applied to the shell has a strong effect on the fragmentation process. 
In cases where shells are not pressure-confined, the shells thicken as they expand and 
hydrodynamic fiows perpendicular to the plane of the shell suppress fragmentation at 
short wavelengths. If the shells are pressure-confined internally and externally, so that 
their thickness remains approximately constant during their expansion, the agreement 
with the analytical solution is better. 
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1 INTRODUCTION 

Expanding shells or bubbles are an ubiquitous feature of 
the interstellar medium (ISM) of the Milky Way and of 
other galaxies. Radio survey s of the Northern ([Heileslll97' ' 
1 19841 ) and Southern skies (jMcClure- Griffiths et alTlm 
have detec ted many t ens of shell-like structures, in HI emis- 
sion, and E hlerova P alous (2005) have used an automated 
search algorithm to find more than 600 in the outer Milky 
Way within the Leiden-Dwingeloo HI survey. They calcu- 
late that the three-dimensional filling factor (defined as 
the volume fraction of the Galaxy occupied) of such shells 
must be at least 5%. HI shells tend to be large (10-1000 
pc) and correspondingl y ol d (1-10 Myr) objects. R ecently 
IChurchwell et al.1 (|200 6h and' Churchwell et al.1 (|200/t ) report 
the detection of ~ 600 infrared shells using the GLIMPSE 
survey. Infrared shells are revealed by dust and PAH emis- 
sion aii^_axe_sm^^ pc) and younger (< IMyr). 

IChurchwell et all (|2006l ) find that the scale height of 
their sample of bubbles is ^ 60pc, less than that of super- 
nova remnants and much less than that of planetary nebu- 
lae, but very similar to the scale height of HII regions. They 
show that the bubbles are likely connected to very young OB 
stars. They find that ~ 25% of their bubbles contain known 
HII regions, and ^ 13% contain open clusters, whereas they 
find only a few to be connected with known supernova rem- 
nants and none with known planetary nebulae. They suggest 
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that the ~ 75% of bubbles not associated with HII regions 
may be powered by stars later than B4 which do not have 
strong enough ionizing fluxes to produce detectable HII re- 
gions, and conclude that the majority of infra-red bubbles 
are driven by feedback from OB stars. 

Such shells have a profound effect on the dynamics of 
the ISM and, in particular, may be respo nsible for trigger- 
ing star fo rmation. N umerous aut hors (e.g. ^Deharveng e t al.l 
2003; Zava gno et al . 2006; Dehar veng et al . 2006, 2008) find 
evidence for populations of young (and oft en massive) stars 
in the peripheries of bubbles, and Ichurchwell et all (|2007[ ) 
find that ^ 18% of their sample of 269 bubbles in the inner 
Milky Way show evidence of triggered star formation. How- 
ever, Church well et al. (2007) mention only the triggering 
mechanism of the bubble overrunning pre-existing clumps 
in the ISM and causing them to collapse. They do not con- 
sider the possibility of the shells themselves becoming gravi- 
tationally unstable and collapsing, even though they observe 
that bubbles often have filamentary structure. 

Numerous authors have studied the gravitational 
fragmentation of expandin g shells using the thin- 
shell approximation, (e.g. Vishniad 1983 : ElmegreenI 

'1994; Whitwo rth et all Il994al lbl: IWiinsch PalousI I2OOII : 
Parmentier 2003), and we will not give a detailed disqui- 
sition here. The analysis consists of studying the stability of 
perturbations of a given angular size in the surface density 
of an infinitesimally-thin shell and hence deriving a disper- 
sion relation, i.e. an expression for the rate of growth uj of 
a perturbation as a function of its angular size 77. It is cus- 
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tomary to denote the angular size of perturbations by 77, but 
we break with this convention and use / instead since our 
analysis consists of decomposing the surface density pertur- 
bations into spherical harmonics, where each spherical har- 
monic mode is characterised by the coefficients / and m {rj 
and / may be used interchangeably for describing structures 
of a given angular size). 

The thin-shell dispersion relation has the form 



2R 



R 



(1) 



where V, R and E are respectively the instantaneous 
velocity, radius and surface mass density of the shell, Cs is 
the sound speed within the shell. This dispersion relation 
applies in the case of a shell expanding in a vacuum and is 
slightly different from that appr opriate for a shel l expanding 
in a constant-density medium (|Elmegreenlll994l ). 

Since, u is the perturbation growth rate, a positive 
value of 00 implies a growing perturbation, a purely real 
negative value implies dissolution and a complex value 
indicates the damped oscillation of a perturbation. 

Despite the heavy reliance of theoretical astrophysics 
on numerical simulations, papers presenting compar- 
isons of different numerical schemes, or even of different 
codes using; the sam e scheme are s till re l atively rare (e.g;. 
Durisen et all Il986l: IPavies et all Il993l: [Bate BurkertI 
19971: lO'Shea et all l2005l: Ide Val-Borro et all I2OO6I : 
Commercon et al.1 I2OO8I : iKitsionas et all I2OO8I ). In this 
paper, we compare not only two different codes, but two 
radically different numerical schemes, namely the Eulerian 
Adaptive Mesh R efinement (AMR) ra ethod (represented by 
the FLASH code. lFrvxell et al.1 ('2OOO')), and the Lagrangian 
Smoothed Particle Hydrodynamics (SPH) technique (rep- 
resented by a variant of the code described in Benzl (|199Q )V 
The gravitational fragmentation of a thin shell is an ideal 
test of the two codes since its relative simplicity allows 
the derivation of a detailed analytic solution. Despite this 
simplicity, it is still a problem of considerable astrophysical 
significance, particularly in the field of triggered star 
formation. This allows us to simultaneously test the two 
numerical methods and to test the results against the 
predictions of the thin-shell approximation and thus to 
explore its relevance to realistic problems. 

The outline of the paper is as follows: in section 2, we 
describe our numerical methods. In section 3, we discuss in 
detail and justify our choice of initial conditions. Section 
4 contains the results of our comparative study, while in 
Section 5 we discuss the astrophysical implications of this 
work. We draw our conclusions in Section 6. 



2 NUMERICAL METHODS 

The purpose of this paper is twofold. Our primary aim is to 
test the reliability of the thin-shell approximation in real- 
istic astrophysical situations. Our methods involve sophisti- 
cated numerical hydrodynamics codes and we make use of 
two very different computational schemes, largely to ensure 
that our results are robust, but also with the secondary aim 
of testing the two codes against each other. 



There are essentially two ways of representing a fluid nu- 
merically. Eulerian schemes impose a stationary mesh on the 
simulation volume which can be refined at will to ensure ad- 
equate resolution of high-density features or regions where 
there are strong gradients in physical quantities, or made 
coarser in regions where such gradients are shallow - such a 
technique is termed Adaptive Mesh Refinement. Conversely, 
the Lagrangian Smoothed Particle Hydrodynamic method 
represents the fluid by discrete objects which can move any- 
where under the influences of the forces acting upon them, 
and whose dimensions are automatically adjusted to give 
the highest resolution to the densest regions of the gas. 

We ch oose to use the pu blicly available FLASH code 
version 2.5 (jFrvxell et aLlbOOOI V It uses the block-structured 
AMR technique imple mented in the PARAMESH library 
(jMacNeice et aLlbood V Its hydrodynamic module uses the 
piecewise parabolic method (jColella & Woodward 1983) 
which incorporates a Riemann solver to compute fluxes be- 
tween individual cells. The code is parallelized using the MPI 
librarjQ. Self gravity is included using a fast tree algorithm 
to be described in a forthcoming paper. 

Our SPH code is a variant of that described in 



Bend (I1990I ) but more recently updated and described in 
Bate et al. | (|l995l l. The fluid equations are solved using the 



SPH formalism, includin g the stand a rd art ificial viscosity 
prescription described in 'Monaghan' ("1992), with a = 1, 
f3 — 2. Self-gravity is implemented by means of a binary tree 
and very high density region s can be replaced b y point-mass 
sink particles, as described in lBate et al.l (|l995l V although we 
do not make use of this facility in this work, since formation 
of sink particles indicates non-linear behaviour and we are 
concerned with testing a theoretical model close to the lin- 
ear regime. We will go on to study non-linear fragmentation 
and to study in detail the mass spectrum produced in later 
papers. 

To ensure a faithful comparison between the two codes, 
we analyse their output using a common technique. We are 
studying the formation of structures in expanding shells, 
which is essentially a two-dimensional problem. We there- 
fore project our shells onto the surface of a sphere at 
a sequence of timesteps in our simulations. These two- 
dimensional surface density maps are then decomposed into 
spherical harmonics using the publicly-available nf ft Fast 
Fourier Transform librarjQ. 

A function / on the surface of a sphere may be repre- 
sented by 



00 m=-\-l 



(2) 



where the Y{^{0^ 0) are spherical harmonic functions, the aj^ 
are the associated coefficients and the angular wavenumber 
/ is given by 



2nR{t) 
A ' 



(3) 



^ The MPI library can be downloaded from 
|http : //www . mcs . anl . gov/research/pro j ect s/mpi/ where 
papers discussing MPI and its applications may also be found. 
^ The nfft library can be downloaded from 
[http: //www-user .tu- Chemnitz .de/^ pott s/nf ft/ 1 
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where A is the wavelength of a given structure. It is cus- 
tomary to work in terms of the angular power spectrum Ci , 
defined by 



Ci 



2/ + 1 



m=-\-l 

E 



(4) 



but we are interested in the surface density perturbations 
themselves, which are described by >/CT- We numerically dif- 
ferentiate the timeseries of a/C^ produced by the nf f t library 
to determine the growth rate of structure as a function of 
/. We can then directly compare our numerically-generated 
dispersion relation a;num(0 — d\n{y^i)/dt with the analytic 
^ani(0 given in Equation [T] 



3 SIMULATIONS 

3.1 Avoiding instabilities 

The primary aim of this paper is to study the gravitational 
fragmentation of expanding spherical shells. These shells 
are usually envisaged as being driven internally by an 
expanding HII region, wind bubble or supernova bubble, 
and expanding into an external medium. The internal 
driving and the external medium introduce two major 
instabilities, respectively the Rayleigh-Taylor and Vishniac 
instabilities, which may complicate the evolution of the 
shell and which we have been careful to avoid. We have 
done this partly to ensure that our comparison with the 
thin-shell approximation is pure and involves only the 
gravitational instability, and partly because of the extreme 
difficulty of simulating such instabilities on a global scale. 

Internal driving of the shell may introduce the 
Rayleigh-Taylor instability, which results from a low- 
density fluid accelerating a high-density fluid. This leads to 
mixing at the interface and the descent of 'fingers' of dense 
fluid into the rarefied fluid. Transverse flows between the 
two fluids then rapidly develop and the Kelvin-Helmholtz 
instability then also appears at the interface. The flow then 
becomes extremely complex and difficult to simulate. To 
avoid the development of the Rayleigh-Taylor instability, 
the shell must be allowed to evolve in freefall with equal 
pressures on the inner and outer surfaces of the shell. 

The sweeping up of an extern al medium can l ead to 
a second instability described by IVishnia3 (|l983l ). The 
Vishniac instability results from the interplay between the 
ram-pressure of the material being collected by the shell, 
which acts on the outer face of the shell and is always 
directed radially inwards towards the shell centre, and 
thermal pressure acting on the inner face of the shell, 
which is always locally perpendicular to the shell surface. 
The different geometries of these pressures drive transverse 
flows in the shell and may accelerate fragmentation. To 
eliminate the Vishniac instability, the ram pressure of the 
external medium must be negligible in comparison to the 
the thermal pressure within the shell itself. 

This study focuses on shells of fixed mass expanding in 
a vacuum, or into a medium of such low density that the 
ram pressure exerted and the mass collected are negligible. 



3.2 Radial density profile of the shell 

We adopt a shell radial density profile such that the shell is 
in instant aneous h y drost atic equilibrium at any given radius. 
Following ISpitzed (|l942l l. the profile in the z-direction of 
a self-gravitating layer extending to infinity in the x- and 
y- directions and centred at 2; = may be derived by the 
equation of hydrostatic equilibrium 

■ (5) 
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dz ^ dz 



with the Poisson equation 
These equations yield 



p{z) = posech 
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where po is the density at ^ = 0. To determine this density 
(and thus the density everywhere), we can integrate through 
the slab to find the surface density: 



tanh 
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Letting zo ^ 00 yields 



po 



2ci 



(8) 



(9) 



The above analysis is valid in plane-parallel geometry. In 
spherical geometry with r being the radial coordinate and 
R being the radius of the densest part of the shell, if the 
thickness of the shell is small compared with its radius, we 
may set z = r — R everywhere in the shell and use Equa- 
tion [7] to calculate the radial density profile of the shell if the 
shell has a finite thickness. The shells described by Equation 
[71 are formally infinitely thick, since the density only falls to 
zero at infinity. However, we consider cold shells which are 
bounded inside and outside by hot rarefied gas exerting a 
finite (although in some simulations, negligibly small) pres- 
sure on the shell. We therefore truncate the density profile 
of the cold gas, giving the shell a finite initial half-thickness 
Zo , such that the cold gas and hot gas are in pressure equilib- 
rium at 2; = =bZo. The relation between the external pressure 
Pext and the value of Zo follow from Equation [71 



p(Zo)c, 



: poCgSech^ 



Zo 



(10) 



recalling that Cs denotes the sound speed in the cold gas 
that makes up the shell. Rearranging Equation [HI gives 
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Inserting Equation \TT\ into Equation \T0\ and using 
sech^ [tanh~-^(x)] = 1 — yields 



2c1po 



and finally 
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Figure 1. Evolution with time of the shell thickness Zq for four 
values of the external confining pressure (solid lines). For compar- 
ison, we also plot the evolution of the shell half-thickness (taken 
to be the thickness at which the density p{z) falls to half its max- 
imum value) for the Pext = 1 x 10"-*^^ dyne cm~^ case (dashed 
line). 



This expression for pQ may be inserted into Equation 1111 to 
produce an expression for the shell thickness in terms of the 
external pressure and the surface density: 



ttG [2Pext + ttGSIJ 



tanh 



2Pext + TtGEI^ 



(14) 



Since the shell expands and then contracts, Szo(^) = 
M/(47ri?(t)2). In Figure [1] we plot the variation with time 
of Zq for a shell with the parameters given above and four 
different values of the external pressure. In the high-pressure 
case (Pext = 10~^^dyne cm~^), the shell initially becomes 
thinner as it expands, the thickness reaching a minimum 
value of ^ O.OSpc. This shell would be very difficult to re- 
solve, requiring a prohibitively large number of SPH parti- 
cles or many levels of mesh refinement. Conversely, in the 
case of a confining pressure Pext = 3 x 10~^^dyne cm~^, 
the shell thickens considerably during its evolution, particu- 
larly in the first 5Myr. We therefore choose an intermediate 
value of the the external pressure, Pext = 10 dyne cm~^ 
for which the shell is sufficiently thick to be resolvable and 
where the thickness of the shell is approximately constant 
for the duration of its evolution. 



3.3 Our initial conditions 

In order to avoid the Rayleigh-Taylor and Vishniac insta- 
bilities, we study the momentum-conserving expansion of a 
shell of fixed mass 2 x IO^Mq whose initial radius is lOpc 
and whose initial outward radial velocity is 2.2km s~^. 



chosen so that the shell will expand to radius of ^ 23 pc 
in a time of ~ 15Myr, before beginning to contract again. 
This velocity is rather low in comparison to those observed 
in real shells, but was selected to ensure that the kinetic 
energy associated with the expansion was less than the 
shell's gravitational binding energy, and to keep the ratio of 
the shell's maximum radius to its initial radius small so that 
adequate numerical resolution could be maintained during 
the whole duration of the shell's evolution. We treat the gas 
in the shell isothermally and set its temperature to lOK. 
The shell is purely momentum-driven, being powered only 
by its initial outward radial velocity. We have calculated 
six models which differ in the external pressure applied to 
the shell and in the perturbations of density and transverse 
velocity imposed. In all cases, the pressures within the void 
inside the shell and in the space outside the shell are equal, 
so that there is no acceleration of the shell due to these 
pressures. 

We consider two cases which we refer to as non- 
pressure-confined and pressure-confined. In an AMR 
simulation, it is not possible to simulate regions of true 
vacuum, so in all of our FLASH simulations, the voids inside 
and outside the shell are in fact filled with very low density 
gas. In the non-pressure-confined AMR simulations, the 
thermal pressure, Pext = 10""^^ dyne cm~^, exerted on the 
shell by this gas is very low and never becomes dynamically 
important - we neglect it entirely in the corresponding SPH 
simulations. In the pressure-confined AMR simulations, the 
temperature of the low-density gas inside and outside the 
shell is higher, so that the Pext = 10 dyne cm~^, which 
ensures that the thickness of the shell stays roughly con- 
stant as it expands. In the corresponding SPH calculations, 
an 'external' pressure term is subtracted from the pressure 
of each SPH particle when calculating pressure gradients 
(so that particles whose internal pressure is less than the 
external pressure effectively exert a negative pressure on 
their neighbours). In all cases (non-pressure-confined and 
pressure-confined), the pressures inside and outside the 
shell are kept equal and constant, so that the shell feels 
no acceleration, and the density of the hot gas used in the 
AMR calculations is so low that the ram pressure it exerts 
is neglible in comparison with its thermal pressure. 

In order to study the growth of particular modes, we 
introduce monochromatic spherical harmonic perturbations 
Spert(^,0) with an amplitude of 10% of the mean shell 
surfa ce density, and a correl ated velocity field derived 
from IWiinsch Paloul (|200lh . We chose to study the 
(/,m) = (12,6) and (35,17) modes. The relative amplitude 
of the surface density perturbation was chosen to be 0.1 to 
ensure that significant fragmentation occurs before the shell 
reaches its maximum radius and begins to contract. The 
thin-shell analysis predicts that these modes should grow. 
We also study shells seeded with white noise and compare 
the dispersion relations generated to that predicted by 
the thin-shell analysis in this case. The parameters of our 
simulations are given in Table [1] 

In SPH simulations of the shells with monochromatic 
perturbations, it was necessary to reduce the noise intrinsic 
to the particle distribution, since the level of the noise in our 
initial conditions was comparable to the amplitude of the 
perturbations we wished to insert. If the SPH particles are 
randomly distributed in the shell, the intrinsic noise level 
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Perturbation type 


Si/So 


VI 




Pext 


/ = 12, m = 6 


0.1 


66 ms~ 


-1 


10- ^'''(0) dynecm-2 


; = 35, m = 17 


0.1 


32 ms- 


-1 


lO-i'^(O) dynecm-2 


white noise / Maxwell velocity 


0.1 


200 ms- 


-1 


10-1^(0) dynecm-2 


/ = 12, m = 6 


0.1 


66 ms~ 


-1 


10-13 dyne cm 


/ = 35, m = 17 


0.1 


32 ms- 


-1 


10-13 dynecm-^ 


white noise / Maxwell velocity 


0.1 


200 ms" 


-1 


10-13 dynecm-^ 



Table 1. Initial conditions for our simulations. The first column details the form of the perturbations in the shell, either a monochromatic 
spherical harmonic perturbation or white noise. The secod column gives the relative amplitude of the surface density perturbations and 
the third gives the amplitude (for the monochromatic runs) or dispersion (for the white noise runs) of the velocity perturbations. The 
fourth column gives the pressure applied to the shell (zero for the first three SPH runs). 



in the surface density is ^ a/JV /A/", where N is the average 
number of particles through which a line of sight extending 
to infinity from the centre of the shell must pass. This may 
be estimated as AT « Npart7T(2{h))'^ / (AttR^) , where ATpart is 
the total number of particles in the simulation and (h) is 
the mean smoothing length (note that SPH particles have 
a radius of 2h). The value of (h) depends on the number 
of particles and how the particles are distributed radially, 
i.e. on the thickness of the shell and on its density profile. 
Using 1.2 X 10^ particles to simulate shells of the radii 
and density profiles described in the preceding section, 
we find that (h) ^ O.lpc, so that ^/lN)/N ^ 0.09. Since 
this is very similar to the amplitude of the monochromatic 
perturbations we wished to introduce in the shell, we took 
steps to reduce the noise level. 

Distributing the particles on a uniform grid (e.g. a 
hexagonal close-packed array) can lead to the appearance 
of grid artefacts. This is not usually the case when the shell 
is sweeping up additional SPH particles, since the swept-up 
particles tend to disrupt the uniform grid in the shell and 
wash out such artefacts, but it is a problem in the simulation 
of our constant-mass shells. We therefore used an iterative 
procedure to reduce the noise level in a shell composed of 
randomly-distributed particles by adjusting the particle 
masses. For each particle z, we drew an imaginary line from 
the shell centre to infinity through the centre of the particle 
and generated a list of all the other particles through which 
this line passed. By integrating through the smoothing 
kernels of these particles (including particle z), we obtain 
the real surface density of the shell at the angular position 
of the particle i, ^l^g^i{Oi, (t)i). We then compare this to the 
desired perturbed surface density at the position of particle 
z, Spert(^^5 0O and multiply the mass of particle i by the 
ratio Spert/Sreai- ordcr to avoid numerical difficulties 
arising from having a large range of particle masses, we 
impose the further constraint that no particle's mass may 
differ from the mean particle mass by more than a factor 
of three, so that the total range of particle masses is less 
than ten. Since each particle contributes to ^ N lines of 
sight, this procedure must be repeated iteratively for all 
particles until the discrepancy between the real and desired 
surface densities is everywhere within some tolerance (we 
chose 0.94 < E^ert/^Lai < 1-06). Note that it was not 
necessary to take such steps in the AMR simulations of 
monochromatically-perturbed shells, since the amplitude 
of the intrinsic noise in the AMR density field is very much 
lower than the 10% perturbation amplitude. 

In order to study fragmentation at all angular scales 



simultaneously, we also performed simulations in which 
white noise was introduced in the initial conditions. We 
constructed these by retaining the intrinsic noise in the SPH 
particle distribution and introducing a three-dimensional 
Maxwellian random velocity field everywhere in the shell, 
with the peak of the velocity distribution occurring at the 
sound speed, Cg. Since Cs « V for most of the shell's 
evolution, this causes no gross distortions of the shell 
surface. The presence of this noise in our simulations can 
be justified on physical grounds by noting that no real 
expanding shell will be uniform and, even if it were, the 
instabilities which we have so carefully avoided in these 
calculations would rapidly introduce inhomogeneities. 

To ensure a legitimate comparison between the simula- 
tions of the shells with noisy initial conditions with the two 
codes, we interpolated the SPH density and velocity fields 
onto a grid which was then used as the initial conditions for 
the corresponding FLASH simulations. We used an SPH 
particle number of 1.2 x 10^ and grid resolution of 640^ 
at the highest refinement level (which contains the whole 
shell) to give the codes the same linear resultion of ^ 0.1 pc 

Our three-dimensional simulations are performed in 
a Cartesian coordinate system in which the centre of the 
shell is at {x,y,z) = (0,0,0) and the 2;-coordinate reverts 
to its usual meaning and should not be confused with the 
2;-coordinate defined in Section 3.2. 



4 RESULTS 

We study the evolution of shells which sweep up no mass, 
but with two very different boundary conditions. We first 
describe simulations in which there is negligible pressure 
inside or outside the shell to confine it. As the shell ex- 
pands (i.e. as its radius increases), it also becomes thicker. 
This is partly due to the decrease in surface density and 
self-gravity due to the geometrical dilution of the shell, 
and partly due to the pressure gradients driving material 
away from the dense core of the shell. We find that this 
thickening has an important influence on the fragmentation 
of the shell in that it suppresses fragmentation at shorter 
wavelengths. This represents a restriction on the validity 
of the thin-shell approximation. To determine if agreement 
with the thin-shell analysis can be achieved by preventing 
the shell thickening, we perform a second set of simulations 
in which we apply a much stronger pressure to the inner 
and outer faces of the shell, such that its thickness remains 
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Figure 2. Hammer projections of the initial surface density perturbations in the simulation of the evolution of the monochromatic 
perturbations (left panel: / = 12, m = 6, centre panel: I = 35, m = 17) and the simulation with random-noise initial conditions (right 
panel). The colour scale represents fractional deviation of the surface density E from the mean value (E). 



approximately constant during its evolution. 



4.1 Non— pressure— confined shells 

The / = 12, m = 6 perturbation which we inserted into 
the shell's surface density is shown in the left panel of 
Figure [2] Its evolution is depicted in Figure [3] and we show 
a slice through the z — plane of the shell at a time of 
8Myr in the left half of Figure 2] showing the gas densities 
and mass fluxes in the AMR and SPH calculations. The 
evolution of this mode is similar to, but slower than, that 
predicted by the thin-shell approximation. In Figures [3{a) 
and [Sljb), we compare the instantaneous and integrated 
analytical and numerical growth rates of the 1 — 12 mode. 
As shown in Figure [S^c), after ~ 5Myr the surface density 
perturbation Si = Si (0)exp(/f), where the fragmentation 
integral /f = a;(t )dt , rises to the mean surface density 
So and later the growth of the mode becomes nonlinear. 
The small oscillation in the surface density perturbation 
in the SPH calculation at t ~ 0.5Myr (also visible in the 
other SPH simulations of monochromatic perturbations) is 
a result of fluctuations in high-frequency modes as a result 
of noise in the surface density. The iterative procedure 
described above cannot completely remove the noise from 
the SPH calculations without producing an unacceptably 
large range of particle masses. However, the high degree of 
consonance between the SPH calculations and the AMR 
runs, which have a much lower level of intrinsic noise, 
implies that the noise in the SPH initial conditions has no 
significant impact on our results. 

The left half of Figure |4] shows that strong radial 
flows cause the shell to thicken, but there are also clear 
azimuthal flows feeding material into the density peaks of 
the / = 12 perturbation, allowing it to grow. We see that 
the density within the fragment is slightly higher in the 
SPH calculation. 

The behaviour of the / = 35 perturbation, whose 
initial form is shown in the centre panel of Figure [2] is 
very different. As shown in Figures [5ja) and[5jb), although 
the / — 35 mode does grow for the first ^ 3Myr of the 
calculation, its amplitude then levels off and begins to 
decrease, dropping below its initial amplitude after ~ 7Myr. 

Slices through the z — plane of the shell, shown in 



the right half of Figure 2] again show flows perpendicular 
to the shell, as seen in the simulations of the / = 12 mode, 
and these flows, together with weak tangential flows, are 
transporting material out of the perturbation, causing it to 
dissolve. 

The right panel of Figure [2] shows the initial surface 
density perturbations in the shell with random noise in 
the initial density fleld. Hammer projections of the surface 
density perturbations and dispersion relations at three 
epochs from the SPH and FLASH calculations are shown 
in Figure (6] The codes agree very well with each other. 
Many features in the surface density perturbation are 
common to the output from both codes and there is also 
good agreement in the forms of the numerical dispersion 
relations, although there is some discrepancy in the growth 
rate of long- wavelength modes at early times. We observe 
that fllamentary structures form on the surface of the 
shell which further fragment into strings of approximately 
circular objects. In both calculations, fragmentation at 
all wavelengths is slower than predicted by the thin-shell 
model and fragmentation at / > 20 is suppressed entirely. 
Note that the sharper peaks at higher /-values are higher 
harmonics due to the non-Gaussianity of the forming frag- 
ments at later times, and do not reflect structure forming 
on these scales. These results are in accord with those of 
our calculations involving monochromatic perturbations, 
which showed that low-frequency modes grow, although 
slower than expected, but that the high-frequency modes 
(/ ^ 20) do not grow at all. 



4.2 Pressure— confined shells 

The thickening of the shell caused by its expansion due 
to the overpressure of the shell gas expanding into the 
surrounding vacuum and by the decreasing surface density 
seems to influence the progress of fragmentation. In our 
second set of simulations, we study the evolution of shells 
perturbed in the same way as in the previous section, but 
pressure-conflned such that the shell thickness remains 
approximately constant during the evolution. 

The evolution of the pressure-conflned shells with the 
/ = 12,m = 6 perturbations is very similar to that of the 
same simulations without pressure conflnement. The AMR 
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(a) Variation with time of uji2 in the SPH 
(sohd blue hne) and AMR (dashed blue 
line) simulations, compared with the an- 
alytical value of uji2 (red line). 
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(b) Variation with time of the integrated 
analytical value of uJi2{t) as derived from 
Equation ^ (magenta line) and the inte- 
grated numerical value derived from the 
SPH (solid blue line) and AMR (dashed 
blue line) simulations. 



Time (Myr) 

(c) Variation with time of the mean sur- 
face density Eq (cyan line), and the per- 
turbed surface density as computed by 
Ei(0)exp(/f) (blue lines) and by (Umax — 
^min)/2 (green lines) from the SPH 
(solid lines) and AMR (dashed lines) sim- 
ulations for the I = 12 perturbation. 



Figure 3. Evolution of the / = 12 perturbation in the non-pressure-confined shell. 
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Figure 4. Slices through the z = plane of the non-pressure-confined shells perturbed with the / = 12 (left) and / = 35 (right) modes 
at a time of 8Myr. Top panels are from the AMR calculations, bottom panels are from the SPH runs. Colours represent gas number 
density, arrows are mass fluxes. 



and SPH simulations are in good agreement and once again 
show that the growth of the mode is slower than predicted 
by the thin-shell analysis. The evolution of the surface 
density perturbation (Figures Efc)) is also very similar 
between the pressure-confined and non-pressure-confined 
simulations. Slices through a section of the z = plane 
show that the confining pressure largely prevents the 



thickening of the shell in the radial direction seen in the 
earlier simulations; azimuthal flows are visible at all stages 
of the simulations. In the left half of Figure [HI we observe 
that, instead of the radial flows transporting material 
away from the middle of the shell, these flows are directed 
towards the middle of the shell, increasing the peak density 
and helping to feed material into the density perturbations. 
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(a) Variation with time of the cjss spher- 
ical harmonic multipole coefficient in the 
SPH (sohd blue line) and AMR (dashed 
blue line) simulations, compared with the 
analytical value of c^ss (red line). 
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(b) Variation with time of the integrated 
analytical value of uo2>b{t) as derived from 
Equation ^ (magenta line) and the inte- 
grated numerical value derived from the 
SPH (solid blue line) and AMR (dashed 
blue line) simulations. 
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(c) Variation with time of the mean sur- 
face density Eq (cyan line), and the per- 
turbed surface density as computed by 
Ei(0)exp(/f) (blue lines) and by (Umax — 
^min)/2 (green lines) from the SPH 
(solid lines) and AMR (dashed lines) sim- 
ulations for the I = 35 perturbation. 



Figure 5. Evolution of the / = 35 perturbation in the non-pressure-confined shell. 




Figure 6. Hammer plots from the SPH and AMR simulations of the non-pressure-confined shell seeded with a Maxwellian random 
velocity field at three epochs (a-c) and comparisons of uj(l) derived from the SPH calculation (red solid line) and from the AMR 
calculation (red dotted line) with the analytical dispersion relation (solid blue line) at the same epochs (d-f). Only the real part of uj is 
shown. 
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(a) Variation with time of uji2 in the SPH 
(sohd blue hne) and AMR (dashed blue 
line) simulations, compared with the an- 
alytical value of uji2 (red line). 




2 4 6 

Time (Myr) 

(b) Variation with time of the integrated 
analytical value of uJi2{t) as derived from 
Equation ^ (magenta line) and the inte- 
grated numerical value derived from the 
SPH (solid blue line) and AMR (dashed 
blue line) simulations. 
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(c) Variation with time of the mean sur- 
face density Eq (cyan line), and the per- 
turbed surface density as computed by 
Ei(0)exp(/f) (blue lines) and by (Umax — 
^min)/2 (green lines) from the SPH 
(solid lines) and AMR (dashed lines) sim- 
ulations for the / = 12 perturbation. 



Figure 7. Evolution of the I = 12 perturbation in the pressure-confined shell. 
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Figure 8. Slices through the z = plane of the pressure-confined shells perturbed with the / = 12 (left) and / = 35 (right) modes at a 
time of 8Myr. Top panels are from the AMR calculations, bottom panels are from the SPH runs. Colours represent gas number density, 
arrows are mass fiuxes. 



This behaviour is due to the constant external pressure 
coming to dominate the decreasing pressure inside the shell, 
as predicted in Figure [1] With the chosen external pressure 
of 10"^^ dyne cm , the shell initially thickens slightly, but 
then decreases in thickness from a time of ^ 2.5Myr. 

By contrast, the / = 35,m = 17 perturbation behaves 
very differently in the pressure-confined shell than in the 



shell with negligible confining pressure. Although initially, 
the growth of this mode is much slower than predicted by 
the thin-shell approximation, the growth rate eventually 
accelerates, at around the time when the shell reaches its 
maximum thickness and begins to get thinner. As with the 
longer-wavelength perturbation, the growth of the / = 35 
mode eventually leads to nonlinear behaviour when the 
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(a) Variation with time of c^ss in the SPH 
(sohd hne) and AMR (dashed hne) sim- 
ulations, compared with the analytical 
value of cjss (red line). 
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(b) Variation with time of the integrated 
analytical value of uo2>b{t) as derived from 
Equation ^ (magenta line) and the inte- 
grated numerical value derived from the 
SPH (solid blue line) and AMR (dashed 
blue line) simulations. 
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(c) Variation with time of the mean sur- 
face density Eq (cyan line), and the per- 
turbed surface density as computed by 
Ei(0)exp(/f) (blue lines) and by (Umax — 
^min)/2 (green lines) from the SPH 
(solid lines) and AMR (dashed lines) sim- 
ulations for the / = 35 perturbation. 



Figure 9. Evolution of the / = 35 perturbation in the pressure-confined shell. 




Figure 10. Hammer plots from the SPH and AMR simulations of the pressure-confined shell seeded with a Maxwellian random velocity 
field at three epochs (a-c) and comparisons of uj{l) derived from the SPH calculation (red solid line) and from the AMR calculation (red 
dotted line) with the analytical dispersion relation (solid blue line) at the same epochs (d-f). Only the real part of uo is shown. 
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surface density perturbation due to the mode comes to 
exceed the mean shell surface density. 

The reason for the markedly different behaviour of this 
mode in the pressure-confined shell is shown in the right 
half of Figure [SI We see that, contrary to the simulation of 
the non-pressure-confined shell, azimuthal flows and radial 
flows towards the middle of the shell feed material into the 
density peaks of the / = 35 perturbation, allowing them 
to grow. We again observe that the peak density is higher 
in the SPH simulation, which is partly due to the SPH 
shell being slightly thinner than the AMR shell, due to the 
different ways in which the external pressure is treated in 
the two codes. 

The results of the two previous simulations imply 
that fragmentation in the pressure-confined shell proceeds 
in a manner more consonant with the predictions of the 
thin-shell approximation. To see if this impression is 
correct, we turn to the simulations of the shell with the 
white noise in the initial density field and the Maxwellian 
velocity distribution. 

In these simulations, we observe again that the agree- 
ment between the two codes is very good and that there is a 
much closer agreement between the numerical calculations 
and the thin-shell approximation. The numerical dispersion 
relations, shown in the lower part of Figure [10] at early times 
show a similar deficit in fragmentation at shorter wave- 
lengths, but the agreement with the theoretical dispersion 
relation improves and, in later stages of the simulations, 
the simulations show fragmentation at all wavelengths pre- 
dicted by the thin-shell approximation. The growth rates 
observed are also in good agreement with the theoretical 
treatment. As a result, structure in the pressure-confined 
shell grows much more rapidly than in the shell without 
a confining pressure (Figure [6]). The fragmentation shows 
the same general form in that it appears to be filamentary 
in nature. It is also clear from a comparison of Figures [6] 
and [To] that the latter exhibits structure at smaller wave- 
lengths, as expected from the numerical dispersion relations. 



5 DISCUSSION 

It is clear from the results presented in Section 4 that 
the predictions of the thin-shell analysis as to the rates 
at which structures of a given wavelength will grow in a 
self-gravitating shell should be treated with some caution. 
The thin-shell approximation treats the shell as infinitesi- 
mally thin, so that the idea of the boundary conditions on 
the inner and outer surfaces of the shell is meaningless. Of 
course, no real shell is infinitesimally thin and we find that 
the boundary conditions applied to the shell surface affect 
the fragmentation very strongly. In a shell expanding in the 
absence of confining pressure, the radial thickening of the 
shell suppresses fragmentation at short wavelengths. The 
radial flows acting away from the middle of the shell dilute 
the azimuthal flows which would otherwise feed density 
perturbations, so that short- wavelength perturbations 
are unable to grow. We therefore find that the thin-shell 
approximation produces a rather poor description of the 
fragmentation occurring in such a shell. The simulations 
show that shells expanding without pressure confinement 



underproduce fragments at short wavelengths. This result is 
potentially interesting in an astrophysical context, since it 
would imply that such shells will produce a top-heavy mass 
function - either a cluster mass function or a stellar IMF, 
depending on the parameters of the shell, in particular the 
total mass. Other effects, such as gas expulsion can also 
infiuence the clu ster mass function and m ake it top-heavy, 
as discussed in jParmentier et a L ^ ('2008'). In either case, 
massive stars will be over-represented with respect to 
low-mass stars. This in turn makes the idea that star 
formation can self propagate more plausible, since triggered 
star formation must produce massive stars in order for this 
to occur. 

Conversely, we find that an expanding shell which is 
pressure-confined internally and externally in such a way 
that the thickness of the shell remains approximately con- 
stant during the shell's expansion fragments in much better 
agreement with the thin-shell approximation. Because the 
shell is prevented from thickening in the radial direction, 
radial fiows in fact act towards the middle of the shell and 
azimuthal fiows dominate over radial fiows in the neighbour- 
hood of perturbations. Hence, surface density perturbations 
are able to grow at all spatial frequencies predicted by 
the thin-shell approximation at rates in good agreement 
with the theoretical predictions. Fragmentation occurring 
in this way should produce an IM F very similar in form 
to the Salpeter IMF ( Paloud lioOTl ) . Furthermore, since the 
surface of the shell becomes deformed as fragments develop, 
the confining pressure accelerates fragmentation. Near the 
forming fragments, components of the pressure forces can 
act tangentially to the shell, driving more material into 
the fragments. We describe this effect as 'pressure-assisted' 
gravitational fragmentation and will study it in detail in a 
subsequent paper. 

The thin-shell analysis assumes that all spatial modes 
evolve independently of each other, so that there is no 
communication or transfer of power between modes. It 
is not obvious that this should be true. We are able to 
achieve good agreement with the thin-shell analysis if the 
thickness of the shell is kept constant during the expansion 
(Figure [TO]) . If mode-mode interaction implies the growth 
of some modes at the expense of others, there should be a 
corresponding change in the form of the dispersion relation. 
Since we do not observe this, we infer that there is no 
interaction of this type in our simulations. 

We have deliberately tried to eliminate both the 
Rayleigh-Taylor and Vishniac instabilities in our calcula- 
tions by studying a shell expansion with no ram-pressure 
and by ensuring that there was no pressure differential 
across the shell. In reality, both these instabilities are likely 
to act on the shell and both may induce it to fragment in a 
manner very different to that produced by the gravitational 
instability acting alone. The classical Rayleigh-Taylor 
instability acts on all spatial lengthscales and it acts fastest 
on the smallest ones, so this instability is likely to generate 
small-scale structure in the shell much faster than the 
gravitational instability is able to. 

In their study of infrar ed shells in the Milky Way, 
IChurchwell et al.1 ()2006l ) and IChurchwell et al.1 (|2007l ) find 
a strong correlation between the radius R of shells and 
their projected thickness AR, such that AR oc R. Since 
shells are seen projected onto the sky, even when spherical 
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symmetry can be assumed, the thickness measured in this 
way is not the same as the three-dimensional thickness 
of the sheh. However, it is very difficult to conceive of an 
evolutionary scheme for the shell in which the apparent 
projected thickness of the shell increases but the true 
three-dimensional thickness does not. This result there- 
fore implies strongly that the real thicknes s es of shells 
increase as they expand. Churchwel l et al.l (|20Q6h and 
IChurchwell et al.l (|200 7) also show that the shells often have 
intricate filamentary substructure. Based on the results 
presented in this paper, we suggest that these shells are 
likely to be undergoing gravitational fragmentation but that 
the thin-shell analysis will not give a reliable description of 
this process, given that the shells are becoming thicker as 
they expand. We have shown that this circumstance leads 
to the suppression of low-mass fragments and will result 
in a top-heavy IMF. The shells observed by Churchwell 
therefore support the notion that star formation can be 
self-propagating, since the production of massive stars is 
essential for this to occur. O bservations of Sh-217 and 
Sh-2 19 (|Deharve ng et al.'"2QQ3 |), RCW 79 (IZavagno et al.1 
I2OO6I ) and Sh-104 ( Deharv eng et al.1 I2OO3I ) reveal that 
expanding HII regions do indeed sweep up shells of dense 
material and trigger the formation of massive objects which 
in turn generate their own ultracompact HII regions and 
may subsequently drive further star formation. However, if 
such shells produce a top-heavy IMF, they cannot make a 
great contribution to the overall stellar population, since 
they do not produce an IMF consistent with that observed 
on large scales. 

This work is of particular relevance in the context 
of the self-enrichment scenario for the the formation of 
globular clusters. The first generation of massive stars in the 
core of a proto-globular cluster explode as supernovae and 
sweep the remaining gas into a dense shell. Since globular 
clusters are approximately spherical and possess smooth 
potentials, the shell is likely to be spherical and relatively 
smooth and can thus in principle be well approximated by 
the thin-shell approximation. 

The fate of such a shell - whether or not it is bound 
and whether or not it fragments to form a second generation 
of stars - is of crucial importance for several outstand- 
ing questions in globular cluster formation. It remains 
poorly understood why, in our galaxy, the metallicity 
distribution function of halo field stars extends down to 
[Fe/H] ~ —5, while the most metal-poor globular clusters 
have [Fe/H] ~ —2.5. In addition, many Galactic and some 
Magellanic C loud g l obula r clusters exhibit abundance 
pecuUarities (|Smithl I2OO6I , a nd references the rein) , or 
multiple main-sequences fe.g. IPiotto et al.ll2007f ). Numer- 
ous authors have suggested that these features could be 
explained by two or more rounds of star formation taking 
place within globular clusters, with the metals from each 
generation polluting the stars of the next. Wolf-Rayet, 
OB stars and AGB stars a re likely to be the chief sources 
of pollution (|Smithl I2OO6I ). so an understanding of the 
mass functi ons of later gener a tions of stars is crucial to 
this model. IParmentier et al showed that it was 

indeed possible for globular clusters t o retain su p ernova 
ej ec ta and the swept-up shell, and IParmentieii (|2004 l) 
and iRecchi &: Danzigerl (|2005h fur ther showed that it was 
possible for the shells to fragment. IParmentierl (|2004 ll used 



the thin-shell analysis to show that smaller numbers of 
first-generation supernovae and higher external pressure in 
the hot protogalactic background both made fragmentation 
of the shell more likely, since the former condition gives 
the shell a smaller initial velocity and the latter ensures 
that the shell slows down more rapidl y once it leaves t he 
protocluster cloud. As we have shown. jParmentierl (|2004l Vs 
assumption of high-pressure slowing the shell down will 
also prevent the shell from thickening and so her results 
derived from the thin-shell analysis are likely to be correct. 
However, it may not always be the case that such shells 
encounter high enough pressures to keep them thin, in 
which case they will underproduce low mass stars. 



6 CONCLUSIONS 

We have conducted numerical simulations of the fragmen- 
tation of expanding gaseous shells in order to compare the 
numerical techniques of SPH and AMR and in order to test 
the theoretical predictions of the well-known thin-shell ap- 
proximation. 

We find that the agreement between the SPH code and 
the FLASH AMR code is very good despite the differences 
in the numerical schemes and in the way in which constant- 
pressure boundary conditions are applied in the two codes. 
The numerically-derived dispersion relations generated by 
the codes agree very well. 

We find that the thin-shell approximation should be 
used with some caution. When dealing with the concept of 
an infinitesimally-thin shell, the concept of boundary condi- 
tions has no meaning, but in any finite-thickness shell, the 
behaviour of the shell surface must be constrained in some 
way. We find that the boundary conditions chosen have a 
strong effect on the fragmentation of the shell. In the ab- 
sence of confining pressure, flows perpendicular to the shell 
surface cause the shell to thicken as it expands, diluting sur- 
face density perturbations. These flows preferentially wash 
out smaller perturbations and thus suppress fragmentation 
at short wavelengths. Under these conditions, the thin-shell 
analysis gives a poor description of the fragmentation pro- 
cess and the shell produces a top-heavy fragment mass func- 
tion due to the deficit of small objects. 

Conversely, if a constant pressure is applied to the inner 
and outer surfaces of the shell such that its thickness does 
not change significantly as it expands, we obtain very good 
agreement with the thin-shell approximation in the sense 
that the analytical predictions of which wavelengths are un- 
stable and which stable match the numerical computations, 
although we generally find that fragmentation at all wave- 
lengths is slower than predicted by the thin-shell model. 

We note that fragmentation of the shells perturbed only 
by random noise (either in the density or velocity fields) re- 
sults in filamentary structures. This effect has been observed 
in simulations of this type before (Dale et al. 2007). Al- 
thou gh they do not mention t he co ncept of shell fragment a - 



tnou gn tney do not mention t ne co ncept 01 sneli iragmenta - 
tion. IChurchwell et al.1 (|2006l ) and IChurchwell et al.( (|2007l ) 



note that many of the infrared bubbles in their observations 
exhibit a similar filamentary structure. 
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